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Abstract: Deconvolution algorithms have been used successfully for optimization/restoration/deblurring 
of astronomical images in a variety of wavelengths, especially in the optical band (e.g., for HST). We 
present here an iterative RiCHARDSON-LuCY type method designed for treatment of images obtained 
with the H.E.S.S. array of ground-based gamma-ray telescopes. Its application to shell-type supernova 
remnant images yields refined details relevant for the study of correlations with other wavelengths, and 
hence for interpretation in terms either of hadronic or leptonic origin of the observed VHE gamma-ray 
emission. 
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Introduction 

All images or maps obtained with a telescope are 
a statistical realization of the real source distri- 
bution degraded by the instrument response and 
with a certain noise: a point source appears thus 
with a spread distribution called the point spread 
function (PSF). So, in order to correcte these ef- 
fects, we have to subtract noise and to deconvolve 
the map with the PSF. Several methods are used 
to do this treatment but they are not adapted for 
VHE gamma-ray astronomy in which the noise 
is comparable to source signal and in which we 
have less statistics. In this paper, we adapt the 
RiCHARDSON-LuCY method, based on a maxi- 
mum likelihood (searching for the most probable 
source distribution) to these conditions. We then 
apply this method to study the shell of the syper- 
nova remnant RXJ1713.7-3946. 



propability density relationship: 



P{0\N,) = 



p{Ns\0)p{0) 

p{Ns) 



(1) 



where p{Ns) is the probability of the observed im- 
age Ns{x,y) and p{0) is the probability of the real 
image 0{x, y), i.e. the source distribution, over all 
possible image realizations. The bayes solution is 
found by maximizing the right part of the equation, 
i.e. by searching the source distribution 0{x,y) 
which maximizes the probability (1). 

The term P{Ns) has a constant value and can be 
ignored during the maximization process, while we 
can assume a uniform probability density for the 
prior p{0). So the maximum probability corre- 
sponds to the maximum likelihood over 0{x,y): 



ML{0) =maxp(Ar^|0) 



(2) 



The classical 
Method 



Richardson-Lucy 



The RiCHARDSON-LuCY method algorithm was 
developped from Bayes's theorem [[6, 5]]. Its ap- 
proach consists of constructing the conditionnal 



The probability P{Ns\0) depends on the object 
and the noise. In our case, they follow Poisson dis- 
tributions. The probability p{Ns\0) is then given 
by: 

-r-r iV (a; u)-'V.(x,y) -iV,(x,y) 



x,V 



Ns{x,y)\ 
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where Ns is the convolution beetween the source 
distribution 0{x, y) and the normahzed point- 
spread function Psf{x, y): 



Ns{x, y) = 0{x, y) * PsJ{x, y) 



(4) 



The maximum UkeUhood can be computed by tak- 
ing the derivative of the logarithm: 



d\np{Ns\0) 
do 



= 



(5) 



Using Picard iteration [[3]], we have: 
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Ns{x,y) 



0{x,y) *Psf{x,y) 



*Psf{x,y) 

(6) 

The maximum likelihood is then obtained when all 
the coefficients a"(a;, y) converge to 1. This algo- 
rithm, commonly used in astronomy, intrinsically 
applies the positivity constraint of the solution and 
conserves flux. But this method doesn't work when 
the noise becomes important. Indeed, in this case, 
it creates after a few number of iterations high fre- 
quency structures which are not physical. We thus 
propose to modify this method in order to work 
with gamma-ray maps which present an important 
background. 

Deconvolution of maps with noise 

Commonly in gamma-ray astronomy, a back- 
ground measurement is made for each observation. 
There are in fact several methods in order to do this 
study (ring background, template background...) 
but a noise map Noff{x, y) is always obtained and 
can be sustracted to the Non{x, y) map to compute 
the excess map. In fact, the image data Non is as- 
sumed to be due to the source distribution 0{x,y) 
convolued by the PSF and to a noise Noff{x,y) 
by: 

Non{x, y) = (0(x, y) * Ps/(x, y))A{x, y) 
+Noff{x,y) 



where A{x, y) is the gamma acceptance-ray in the 
feld of view of the telescopes. 

Then we can construct a conditionnal propability 
density relationship, in order to take into account 
the signal and noise distributions: 

p{Non\Non)piNoff\Noff)^ 

p{Non\Non)p{Non) ^ P{Noff\NofMNoff) 



p{No 



p{No 



ff) 



Using the same assumptions than in the classical 
method, the solution is obtained for the maximum 
likelihood over 0{x,y) and over Noff- These 
maximum likelihoods can be computed by taking 
the derivative of the logarithm: 



91n [p{Non\Non) P{Noff\N~ff)] 



do 



= 



^ 91n[p(No„|Non)p(Noff'|N;ffO] „ 
and -= — 

aNoff 

The maximum likelihood over Noff, is obtained 
for: 

Noffix,y) ^ Noff{x,y) (7) 

So Noff can be directly remplaced by Noff{x, y) 
for the maximum likelihood over 0(x, y). We find 
then a new iteration coefficient a"(x, y) in Picard 
iteration given by: 



a(")(2;,2/) 



Nar,{x,v)A{x,y) 



{0'-^){x,y)*Psf(x,y))A(x,y)+N„fs{x,y) 



* Psf{x, y) 



A{x,y) *Psf{x,y) 

In this way, we take into account both noise and 
acceptance while allowing negative values for the 
source distribution. To avoid this problem, we can 
impose 0{x,y) positive in each pixel during the 
Picard steps. This condition brakes however the 
photometry, but we have show on simulations that 
this effect is very weak (< 1%). 



When does the modified Richarsond- 
LuCY method converge ? 

There are several criteria used to test the con- 
vergence of the deconvolution method (bootstrap 
method [[!]], comparison with realistic simula- 
tions...). Here we propose a simple approch where 
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we stop the deconvolution algorithm, pixel-by- 
pixel, when the coefficient a"-{x, y) reaches 1 with 
a choosen limit. This Umit depends on the PSF and 
the bin size. So we have to adjust it to the map that 
we want to analyse. The next paragraph presents 
this work on the supernova remnant: RXJ1713.7- 
3946. 

Application to the supernova remnant: 
RXJ1713.7-3946 

The figure 2 (left) presents the excess map of 
RXJ1713.7-3946 obtained by the H.E.S.S. tele- 
scopes [[2]]. The map binning is of 0.02° and 
the PSF is modelized by simulations through a 2- 
dimensionnal Gaussian of ct ~ 0.1°. 
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Estimation of the optimal iteration number 

The supernova remnant shell projected on the sky 
appears as a fraction of a ring with several hot 
spots. To estimate the optimal iteration number, we 
simulate Gaussian- and ring-shaped sources with 
different sizes (from 0.5 to 2 PSF size) and for sev- 
eral intensity values (yielding a significance be- 
tween 10 to 40cr) and with the same binning as 
data. For these simulations, we seach for the best 
iteration number with a limit on a"(x, y) (see pre- 
vious paragraph) and evaluate the spatial resolution 
by fitting the source size after deconvolution. The 
plots given in 1 presents these resolutions with a 
limit of 1%. 

This method makes converge the deconvolution to 
a limit of resolution depending on the significance 
and the source size. The higher the significance 
or the bigger the source, the better the resolution 
is. These tests show also that the deconvolution is 
a slow process after several iterations. So, if we 
don't apply exactly the optimal iteration number, it 
doesn't modify significantly the deconvolued map. 

The deconvolution of RXJ1713.7-3946 

The RXJ1713.7-3946 map through H.E.S.S. obser- 
vations yields a significance of SOcr. The source 
morphology is close to a fraction of a ring of width 
K, 0.2°. Applying the deconvolution with 120 it- 



Figure 1: Source size error as a function of in- 
tensity for Gaussian- and ring-shaped sources with 
different sizes (from 0.5 to 2 PSF size). 

erations (according to fig.l) we obtain the map 
drown in fig. 2 right. 

The deconvolution improves clearly the corre- 
lation between the TeV emission and the X-ray 
emission (1-3 keV) as measured by ASCA [[4]]. 
Figure 3 presents the improvement by plotting the 
TeV emission as a function of X-rays emission 
from the ASCA map for each bin in the east part 
of the map. These correlations are important to 
constrain models. 

In leptonic models, synchrotron emission of elec- 
trons generates X-ray and inverse Compton with 
photon target (optical star light, CMB, IR) gener- 
ates the TeV emission. We then expect strong cor- 
relations. In the hadronic modeld, Tev emission is 
due to the 7r° decay generated by the interactions 
between accelerated hadrons and nuclei in the site. 
These interaction generate electron which could be 
responsible to the X-ray emission. But, at the mo- 
ment, this process has to be check by modelization. 
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Figure 2: Excess (left) and deconvolved (right) map of RXJ17 13.7-3946 observed by H.E.S.S. The excess 
map is smoothed by a gaussian of 0.04° in order to see the structure of the supernova shield. In both maps, 
the ASCA experiment (1-3 keV) contours [[4]] are overdrawn. 
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Figure 3: Correlation between the TeV emission 
from the H.E.S.S. map and the X-rays (1-3 keV) 
from the ASCA map [[4]]. 

Conclusion and future prospects 

Our method, based on a classical algorithm for 
treatment of images in the optical band, allows us 
to improve the angular resolution of gamma-ray as- 
tronomy by making a statistical correction of the 
instrumental spread. This analysis can be already 
used to study the morphology of extended sources. 
In a futur step, the background treatment should be 
improved during the deconvolution process and the 
limits of the algorithm should be estimated. 
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